/* Gionata Massi <gionata.massi@gmail.com>
 * Dicembre 2006- Dicembre 2009
 *
 * File: optimizer.cpp
 *
 * Questo file fa parte di 'aga'
 */

/*!
  \file		optimizer.cpp
  \brief	Solutore MOSP mediante algoritmo genetico.
  \author	Gionata Massi
  \version	2009-12-15
  \date		2009-12-15

  Il nucleo dell'Algoritmo Evolutivo, la generazione di popolazione iniziale e
  nuove generazioni, l'ordinamento della popolazione in base alla fitness.
 */

#include "aga.h"
#if defined (DLL_EXPORT)
#	define NDEBUG
#endif
#include <assert.h>

#define LOG_FN "stderr"

//! Lista delle euristiche per generare alcune soluzioni iniziali.
/*!
   Contiene l'indirizzo delle funzioni euristiche per la generazione di alcune
   soluzioni iniziali pi&ugrave; un puntatore a NULL per permettere il calcolo
   della lunghezza
 */

HeuristicFunction HeuristicsList[] = {
	InsertionCost,
	//InsertionVehicles,
	CompleteSequence,
	NULL
};


inline void evaluate_population(GeneticData * data);
inline void update_k_lexicography(GeneticData * data);
inline int perform_local_search(GeneticData * data);

/*! Valuta gli indici di performance di ogni individuo.
	Aggiorna gli upper e lower bound degli indici di performance per la
	popolazione della generazione corrente.
	Aggiorna i pesi per la valutazione della fitness.
	Aggiorna la fitness di ogni individuo.
*/
inline void evaluate_population(GeneticData * data)
{

	data->LB[TOTAL_COST] = 0;
	data->LB[VEHICLES_NUMBER] = 0;
	data->UB[TOTAL_COST] = INT_MAX;
	data->UB[VEHICLES_NUMBER] = data->CustomerSet.size();

	// Valutazione indici di performance
	for (int nInd = 0; nInd < data->NumIndividuals; ++nInd) {

		if (!data->People[nInd].performance_indicators_evaluated) {
			evaluate_performance_indicators(data,
			                                data->People[nInd].
			                                Sequence,
			                                &(data->People[nInd].
			                                  Score));
			data->People[nInd].performance_indicators_evaluated = 1;
		}
		for (int p = 0; p < parameters; p++) {
			if (data->LB[p] > data->People[nInd].Score.measure[p]) {
				data->LB[p] =
				    data->People[nInd].Score.measure[p];
				if (p == TOTAL_COST
				    && data->People[nInd].Score.measure[p] <
				    data->TheBestOne.Score.measure[p]) {
					indcpy(data, &(data->TheBestOne),
					       &(data->People[nInd]));
					data->stat.last_best_generation[p]=
					    data->CurrentGeneration;
					data->stat.ls_found_opt = 0;
				}
				if (p == VEHICLES_NUMBER
				    && data->People[nInd].Score.measure[p] <
				    data->BestVehicleInd.Score.measure[p]) {
					indcpy(data, &(data->BestVehicleInd),
					       &(data->People[nInd]));
				}
			}
			if (data->People[nInd].Score.measure[p] > data->UB[p])
				data->UB[p] =
				    data->People[nInd].Score.measure[p];
		}
	}
	update_k(data->k, data->Weight[data->CurrentPhase], data->LB, data->UB);
	for (int nInd = 0; nInd < data->NumIndividuals; ++nInd) {
		fitness_function(data, &(data->People[nInd].Score), data->k,
		                 data->LB, data->UB);
#ifdef DEBUG
		if (data->People[nInd].Score.score < -0.0000001)
			assert(data->People[nInd].Score.score >= -0.0000001);
		if (data->People[nInd].Score.score > 1.0000001)
			assert(data->People[nInd].Score.score <= 1.0000001);
#endif
	}
	/* Ordinamento */
	SortPopulation(data->People, data->NumIndividuals,
	               data->NumIndividuals);
}

inline void update_k_lexicography(GeneticData * data)
{
	/* Imposta i pesi delle componenti della fitness per la generazione
	 * degli individui euristici. */
	data->Weight[PHASES][TOTAL_COST] = 0.999;	// 0);
	data->Weight[PHASES][VEHICLES_NUMBER] = 0.001;	//0);

	update_k(data->k, data->Weight[PHASES], data->LB, data->UB);
}

inline int perform_local_search(GeneticData * data)
{
	int ls_trials = 0;
	int i = 0;
	for (; i < data->LocSearchNumBestInd; ++i) {
		if (!data->People[i].LS_performed[LS_2Opt_INDEX]) {
			int ls_depth;
			ls_depth = ls_2opt_ind(data, i);
			if (data->stat.ls_max_depth < ls_depth)
				data->stat.ls_max_depth = ls_depth;
			ls_trials++;
			/*
			if (data->BestMosInd.Score.measure[MOS] <=
			    data->LOWER_BOUND) {
				goto TERMINATE_LS;
			}
			*/
		}
	}
	for (; i < data->LocSearchNumTotInd; ++i) {
		int random_idx =
		    rand() % ((int)(data->NumIndividuals * data->FirstParent));
		if (!data->People[random_idx].LS_performed[LS_2Opt_INDEX]) {
			ls_2opt_ind(data, random_idx);
			ls_trials++;
			/*
			if (data->TheBestOne.Score.measure[TOTAL_COST] <=
			    data->Lower_Bound) {
				goto TERMINATE_LS;
			}
			*/
		}
	}

TERMINATE_LS:
	return ls_trials;
}

//! Il nucleo dell'algoritmo genetico
/*!
  \param data	Puntatore alla struttura contente i dati dell'algoritmo
		genetico
*/

int GeneticOptimizer(GeneticData * data)
{
#pragma region Dichiarazione variabili automatiche e inizializzazioni
	int nIndividuals = data->NumIndividuals;
	/* Numero di geni */
	int nGenes = data->NumGenes;
	// /* Indice dell'individuo corrente */
	// int nInd;
	/* Fitness del miglior individuo della generazione precedente */
//      Ranking oldBestValue;
	int termination;

	/* La nuova popolazione, che sostituira' la vecchia */
	Population pNew = AllocatePopulation(nIndividuals, nGenes);

	/* Controlla che l'allocazione di memoria sia riuscita */
	if (!(pNew))
		exit(1);
	operators_initializer(data);
#ifdef VERBOSE_DEBUG
	printf("operators_initializer\n");
	fflush(NULL);
#endif

	localsearch_initializer(data);
#ifdef VERBOSE_DEBUG
	printf("localsearch_initializer\n");
	fflush(NULL);
#endif

	data->stat.elaps[WHOLE_ALG].start = clock();
	data->beginTime = time(NULL);
#ifdef LOG
	log_ini(data, LOG_FN);
#endif
#pragma endregion

	data->TheBestOne.Score.measure[TOTAL_COST] = INT_MAX;
	update_k_lexicography(data);

#pragma region Generazione Popolazione Iniziale
	/* Generazione Popolazione Iniziale */
	data->stat.elaps[INITIAL_POPULATION].start = clock();
	InitialPopulation(data);
	data->stat.elaps[INITIAL_POPULATION].end = clock();
#ifdef VERBOSE_DEBUG
	printf("Pop iniz\n");
	fflush(NULL);
#endif
#pragma endregion

	/*
	if (data->TheBestOne.Score.measure[TOTAL_COST] <= data->MOS_Lower_Bound) {
		data->stat.elaps[REPRODUCTION_CYCLE].start =
		    data->stat.elaps[REPRODUCTION_CYCLE].end = clock();
		termination = AGA_EMOSLB;
		goto TERMINATE;
	}
	*/
#pragma region Imposta i pesi per fase 0
	/* Imposta i pesi delle componenti della fitness per la prima fase */
	//if (data->FitFunPtr == FitFunMOSP) {
	/* La prima fase */
	data->CurrentPhase = 0;
	update_k(data->k, data->Weight[data->CurrentPhase], data->LB, data->UB);
	// } else {exit (0123);}
	// update_k(data);
#ifdef VERBOSE_DEBUG
	printf("up k\n");
	fflush(NULL);
#endif
#pragma endregion

	evaluate_population(data);
#ifdef VERBOSE_DEBUG
	printf("evaluation and sort\n");
	fflush(NULL);
#endif

	data->stat.elaps[REPRODUCTION_CYCLE].start = clock();

#pragma region Ciclo evolutivo (riproduttivo+raffinamento)
	for (data->CurrentGeneration = 0;
	     !(termination = stopEvolution(data)); data->CurrentGeneration++) {

#ifdef VERBOSE_DEBUG
		printf("\tGenerazione: %d\n", data->CurrentGeneration);
		fflush(NULL);
#endif

#pragma region Controllo cambio fase e aggiornamento pesi
#if (PHASES > 1 || !defined(PHASES))
		/* Setta i pesi per il calcolo della fitness */
		if (data->CurrentGeneration ==
		    (int)(data->PhaseHop[0] * data->NumGenerations)) {

			for (int p = 0; p < parameters; p++) {
				data->stat.best_found_at_phase[data->CurrentPhase +
				                               2][p] =
				                                   data->TheBestOne.Score.measure[p];
			}

#ifdef LOG
			dumpPopulation(data, 0);
#endif
			/* La seconda fase */
			data->CurrentPhase = 1;
			update_k(data->k, data->Weight[data->CurrentPhase],
			         data->LB, data->UB);
			data->ProbMutation = data->InitialMutation;
			/* for (nInd = 0; nInd < nIndividuals; ++nInd)  // gli indici di performance dovrebbero essere validi
			   data->People[nInd].performance_indicators_evaluated = 0; */
			evaluate_population(data);
			if (data->LocSearchFreq > -1)
				perform_local_search(data);
		}
#endif
#if (PHASES > 2 || ! defined(PHASES))
		else if (data->CurrentGeneration ==
		         (int)(data->PhaseHop[1] * data->NumGenerations)) {

			data->stat.best_found_at_phase[data->CurrentPhase +
			                               2][TOTAL_COST] =
			                                   data->TheBestOne.Score.measure[TOTAL_COST];

#ifdef LOG
			dumpPopulation(data, 0);
#endif

			/* La terza fase */
			data->CurrentPhase = 2;

			update_k(data->k, data->Weight[data->CurrentPhase],
			         data->LB, data->UB);
			data->ProbMutation = data->InitialMutation;
			/* for (nInd = 0; nInd < nIndividuals; ++nInd)
			   data->People[nInd].performance_indicators_evaluated = 0; */
			evaluate_population(data);
			if (data->LocSearchFreq > -1)
				perform_local_search(data);
		}
#endif
#pragma endregion

#pragma region Genera la nuova popolazione e rimpiazzo generazionale
		/* Genera la nuova popolazione */
		NewPopulation(data, pNew);
		/* La nuova popolazione sostituisce la precedente */
		{
			Population tmp = data->People;
			data->People = pNew;
			pNew = tmp;
		}
#pragma endregion

		evaluate_population(data);

		/* Termina se la soluzione e' certificata ottima */
		/*
		if (data->TheBestOne.Score.measure[MOS] <=
		    data->MOS_Lower_Bound)
			goto TERMINATE;
		*/

#pragma region Un poco di scolarizzazione
		if ((data->LocSearchFreq > 0)
		    && ((data->CurrentGeneration % data->LocSearchFreq) == 0)
		    && (data->CurrentGeneration != 0)) {
			perform_local_search(data);
			SortPopulation(data->People, data->NumIndividuals,
			               data->NumIndividuals);
		}
#pragma endregion

		/* Aggiorna miglior individuo in struct */
		// data->BestValue = data->People[0].Score;
		bool differ = false;
		for (int i = 0; i < data->NumGenes; i++) {
			if (data->TheBestOne.Sequence[i] !=
			    data->People[0].Sequence[i]) {
				differ = true;
				break;
			}
		}
		if (differ) {
			indcpy(data, &(data->TheBestOne), &(data->People[0]));
			double new_mp = data->ProbMutation * 1.1;
			if (new_mp < MUTATION_THRESHOLD)
				data->ProbMutation = new_mp;
			else
				data->ProbMutation = MUTATION_THRESHOLD;
		}
	}
#pragma endregion

	data->stat.elaps[REPRODUCTION_CYCLE].end = clock();

	for (int p = 0; p < parameters; p++) {
		data->stat.best_found_at_phase[data->CurrentPhase +
		                               2][p] =
		                                   data->TheBestOne.Score.measure[p];
	}

#ifdef LOG
	dumpPopulation(data, 0);
#endif

	if (data->LocSearchFreq > -2) {
		update_k_lexicography(data);
		perform_local_search(data);
	}

TERMINATE:
	data->endTime = time(NULL);
	// Libera un po' di memoria allocata
	FreePopulation(&pNew);
	operators_finalizer();
	localsearch_finalizer();
	// */

	data->stat.elaps[WHOLE_ALG].end = clock();
#ifdef LOG
	log_fini(data);
#endif

	return termination;
}

//! Genera una nuova popolazione e la memorizza in \p pNew.
/*!
  Nel calcolo di nMatch si esegue l'AND con -2 per essere certi che sia pari.
  tmp_1 e tmp_2 sono necessari alla routine di crossover. Dato che tmpVect ha
  nGenes elementi, tmp_1 e tmp_2 sono ottenuti da tmpVect. Cio' non crea
  problemi se 2 * sizeof(unsigned char) <= sizeof(int).

  \param data		Puntatore ai dati dell'algoritmo genetico
  \param[out] pNew	Puntatore alla nuova popolazione
  \param tmpVect
*/
void NewPopulation(GeneticData * data, Population pNew)
{
	/* Indice genitori */
	int idx_parent_1;
	int idx_parent_2;
	/* Genitori */
	Individual *parent_1;
	Individual *parent_2;
	/* Numero di geni */
	int nGenes = data->NumGenes;
	/* Numero di individui */
	int nIndividuals = data->NumIndividuals;
	/* Numero di individui che subiranno una mutazione */
	int nMutation = (int)(data->ProbMutation * nIndividuals);
	/* Numero di individui che genereranno figli, mediante Crossover */
	int nMatch = (int)(data->ProbCrossover * nIndividuals) & (-2);
	/* Contatore degli individui generati */
	int newPeople;
	/* Indice per cicli for */
	int i;

	if (nMutation <= 0)
		nMutation = 1;

	/* Elitismo: porta i migliori individui nella nuova generazione */
	for (i = 0; i < data->Elite; i++) {
		/*memcpy(pNew[i].Sequence, data->People[i].Sequence, nGenes * sizeof(Gene)
		   );
		   for (int j = 0; j < LS_TYPES; j++)
		   pNew[i].LS_performed[j] = data->People[i].LS_performed[j];
		   pNew[i].performance_indicators_evaluated = 1; */
		indcpy(data, &(pNew[i]), &(data->People[i]));
		data->People[i].survival++;
	}

	/* Genera gli individui successivi ai migliori vecchia popolazione */
	for (newPeople = i; newPeople < nIndividuals;) {
		/* Sceglie il primo individuo fra i migliori */
		idx_parent_1 =
		    ChooseParent(data, 0,
		                 (int)(data->FirstParent * nIndividuals));
// (int)(data->FirstPartent * nIndividuals) + 1);

		/* Sceglie il secondo, diverso dal primo, fra tutta la popolazione */
		do {
			idx_parent_2 = ChooseParent(data, 1, nIndividuals);
		} while (idx_parent_1 == idx_parent_2);

		/* Imposta i puntatori ai genitori */
		parent_1 = &(data->People[idx_parent_1]);
		parent_2 = &(data->People[idx_parent_2]);

#ifdef _DEBUG
			for (int i = 0; i < data->NumGenes; i++) {
				assert(data->People[idx_parent_1].Sequence[i] != 0);
				assert(data->People[idx_parent_2].Sequence[i] != 0);
			}
#endif

		/* Finche' non sono stati eggettuati tutti i cx previsti */
		if (newPeople < nMatch) {
			/* Sostituisce nMatch individui con altrettanti figli
			 * provenienti da Crossover */
			Crossover(parent_1, parent_2,
			          &(pNew[newPeople]),
			          &(pNew[newPeople + 1]), data);
			assert(pNew[newPeople].
			       performance_indicators_evaluated == 0);
			assert(pNew[newPeople + 1].
			       performance_indicators_evaluated == 0);
#ifdef _DEBUG
			for (int i = 0; i < data->NumGenes; i++) {
				assert(pNew[newPeople].Sequence[i] != 0);
				assert(pNew[newPeople + 1].Sequence[i] != 0);
			}
#endif
			newPeople += 2;
		} else {
			/* Gli individui selezionati passano nella nuova
			 * popolazione, o se gia' sopravvisuti passano
			 * gli individui piu' prossimi */
			int start_individual = idx_parent_1;
			bool other_individual_to_check =
			    data->People[idx_parent_1].survival;
			while (other_individual_to_check) {
				idx_parent_1++;
				if (idx_parent_1 == data->NumIndividuals) {
					idx_parent_1 = 0;
				}
				other_individual_to_check =
				    (data->People[idx_parent_1].survival
				     && idx_parent_1 != start_individual);
			}

			indcpy(data, &(pNew[newPeople]),
			       &(data->People[idx_parent_1]));
			data->People[idx_parent_1].survival++;
#ifdef _DEBUG
			for (int i = 0; i < data->NumGenes; i++) {
				assert(data->People[idx_parent_1].Sequence[i] != 0);
			}
#endif

			++newPeople;

			if (newPeople < nIndividuals) {
				start_individual = idx_parent_2;
				other_individual_to_check =
				    data->People[idx_parent_2].survival;
				while (other_individual_to_check) {
					idx_parent_2++;
					if (idx_parent_2 ==
					    data->NumIndividuals) {
						idx_parent_2 = 0;
					}
					other_individual_to_check =
					    (data->People[idx_parent_2].survival
					     && idx_parent_2 !=
					     start_individual);
				}
				indcpy(data, &(pNew[newPeople]),
				       &(data->People[idx_parent_2]));
				data->People[idx_parent_2].survival++;
				++newPeople;
#ifdef _DEBUG
			for (int i = 0; i < data->NumGenes; i++) {
				assert(data->People[idx_parent_2].Sequence[i] != 0);
			}
#endif
			}
		}
	}

	/* Mutazione */

	/* Sostituisce nMutation individui con mutazioni di se stessi
	 * cioe' lo scambio di due geni secondo RX o 2-opt */
	for (newPeople = 0; newPeople < nMutation; ++newPeople) {
		/* SCAMBIO Reciprocal eXchange casuale */
		/* Per evitare peggioramenti dei migliori, non si applica ad essi */
		int random_idx =
		    rand() % (nIndividuals - data->Elite) + data->Elite;
		Individual *someone = &(pNew[random_idx]);
		int mutPoint1 = rand() % (nGenes);
		int mutPoint2 = rand() % (nGenes);
		if (mutPoint1 == mutPoint2)
			mutPoint2 = (mutPoint1 + 1) % (nGenes);

		Gene g = someone->Sequence[mutPoint1];

		someone->Sequence[mutPoint1] = someone->Sequence[mutPoint2];
		someone->Sequence[mutPoint2] = g;
		for (int ls = 0; ls < LS_TYPES; ls++)
			someone->LS_performed[ls] = 0;
		someone->performance_indicators_evaluated = 0;
#ifdef _DEBUG
		for (int i = 0; i < data->NumGenes; i++)
			if(someone->Sequence[i] == 0)
				exit(1);
#endif

	}

//      while ((rand() % 100) < 20) {
//              int random_idx = rand() % (nIndividuals / 2);
////            Individual *someone = &(pNew[random_idx]);
////            ls_2opt_sd(someone->Sequence, data);
//              ls_2opt_sd_ind(data, random_idx);
//      }
}

/**
  * Effettua un banale bubble sort della popolazione secondo il ranking.
  *
  * TODO: Sarebbe opportuno un HEAPSORT o almeno uno Shacker SORT
  */
void SortPopulation(Population People, int order, int populationSize)
{
	int swap = 1;

	for (int i = 0; i < order && swap; i++) {
		swap = 0;

		for (int j = populationSize - 2; j >= i; j--) {
			if (People[j].Score.score > People[j + 1].Score.score) {
				Individual tmp_ind;

				tmp_ind = People[j];
				People[j] = People[j + 1];
				People[j + 1] = tmp_ind;

				swap++;
			}
		}
	}
}

//! Crea una popolazione iniziale di nIndividuals
/*!
  Crea una popolazione iniziale di nIndividuals ognuno dei quali ha nGenes
  geni

  \param data	La struttura contente i parametri dell'algoritmo
		genetico
*/
void InitialPopulation(GeneticData * data)
{
	/* Un individuo */
	Individual *onebody;
	/* Numero di individui di cui e' composta la popolazione */
	int nIndividuals = data->NumIndividuals;
	/* Numero di geni */
	int nGenes = data->NumGenes;
	/* Rapporto di popolazione euristica */
	double pH = data->HeuristicsProportion;
	/* Variabili di appoggio o contatori */
	int x, g, temp, next_individual;
	/* Numero di euristiche disponibili */
	int heuristics = 0;
	/* Popolazione di individui da generare con euristiche */
	Population greedyPop = NULL;

	int NumGreedySolutions;
	/* Memoria da allocare per euristiche costruttive */
	unsigned *sequence = NULL;
	// double *score = NULL;
	// int *ordered_sequence = NULL;

	next_individual = 0;

#pragma region Generazione di tutti gli individui da alg euristici costruttivi
	/* Va generata popolazione da euristiche? */
	if (pH > 0.0) {
		data->CurrentPhase = -2;
		/* Dimensione dell'array da passare a euristiche costruttive */
		int dim;

		/* Conta quante funzioni euristiche vengono impiegate */
		for (int i = 0; HeuristicsList[i]; i++)
			heuristics++;

		/* Alloca strutture dati d'appoggio */
		dim = data->NumGenes; //4 * nGenes + 2 * data->NumPieces;
		sequence = (unsigned *)calloc(1, dim * sizeof(unsigned));
		NumGreedySolutions = heuristics * nGenes;
		greedyPop = AllocatePopulation(NumGreedySolutions, nGenes);

		/* Generazione delle sequenze greedy */
		/* Usa ogni euristica disponibile */
		for (g = 0; g < heuristics; g++) {
			/* Seleziona l'euristica */
			HeuristicFunction h = HeuristicsList[g];
			/* A partire da ogni schema */

			for (x = 0; x < nGenes; x++) {
				h(data, sequence,
				  (greedyPop[g * nGenes + x].
				   Sequence), x + 1);
			}
		}
		// Tutte le soluzioni greedy sono state generate
		free(sequence);

		// butta le soluzioni uguali

		int removed = 0;
		for (g = heuristics - 1; g > 0; g--) {
			bool done = false;
			for (x = 0; x < nGenes && !done; x++) {
				bool equal = true;
				for (int k = 0; k < data->NumGenes; k++) {
					if (greedyPop[g * nGenes + x].
					    Sequence[k] !=
					    greedyPop[(g - 1) * nGenes + x +
					              removed].Sequence[k]) {
						equal = false;
						break;
					}
				}
				if (equal) {
					Individual tmpInd =
					    greedyPop[g * nGenes + x];
					for (int i = g * nGenes + x;
					     i < NumGreedySolutions - 1; i++)
						greedyPop[i] = greedyPop[i + 1];
					greedyPop[NumGreedySolutions - 1] =
					    tmpInd;
					NumGreedySolutions--;
					x--;
					removed++;
				}
				if ((g * nGenes + x) == NumGreedySolutions - 1) {
					done = true;
					break;
				}
			}
		}


		int solution_to_copy = (int)
		                       (data->NumIndividuals * data->HeuristicsProportion);
		if (solution_to_copy > NumGreedySolutions)
			solution_to_copy = NumGreedySolutions;
		if (solution_to_copy <= 0)
			solution_to_copy = 1;
		if (solution_to_copy > data->NumIndividuals)
			solution_to_copy = data->NumIndividuals;

		data->LB[TOTAL_COST] = 0;
		data->LB[VEHICLES_NUMBER] = 0;
		data->UB[TOTAL_COST] = INT_MAX;
		data->UB[VEHICLES_NUMBER] = data->CustomerSet.size();


		for (int heu_ind = 0; heu_ind < NumGreedySolutions; heu_ind++) {
			evaluate_performance_indicators(data,
			                                greedyPop[heu_ind].
			                                Sequence,
			                                &(greedyPop[heu_ind].
			                                  Score));
			greedyPop[heu_ind].performance_indicators_evaluated = 1;
			for (int p = 0; p < parameters; p++) {
				if (data->LB[p] > data->People[heu_ind].Score.measure[p]) {
					data->LB[p] =
					    data->People[heu_ind].Score.measure[p];
					if (p == TOTAL_COST
					    && data->People[heu_ind].Score.measure[p] <
					    data->TheBestOne.Score.measure[p]) {
						indcpy(data, &(data->TheBestOne),
						       &(data->People[heu_ind]));
						data->stat.last_best_generation[p]=
						    data->CurrentGeneration;
						data->stat.ls_found_opt = 0;
					}
					if (p == VEHICLES_NUMBER
					    && data->People[heu_ind].Score.measure[p] <
					    data->BestVehicleInd.Score.measure[p]) {
						indcpy(data, &(data->BestVehicleInd),
						       &(data->People[heu_ind]));
					}
				}
				if (data->People[heu_ind].Score.measure[p] > data->UB[p])
					data->UB[p] =
					    data->People[heu_ind].Score.measure[p];
			}
		}
		update_k_lexicography(data);
		for (int heu_ind = 0; heu_ind < NumGreedySolutions; heu_ind++) {
			fitness_function(data, &(greedyPop[heu_ind].Score),
			                 data->k, data->LB, data->UB);
		}

		/* Riodina i primi \p solution_to_copy elementi */
		SortPopulation(greedyPop, solution_to_copy, NumGreedySolutions);

		/* Copia i migliori individui generati con algoritmi costruttivi */

		for (int k = 0;
		     next_individual < solution_to_copy;
		     next_individual++, k++) {

			indcpy(data, &(data->People[next_individual]),
			       &(greedyPop[k]));
			data->People[next_individual].
			performance_indicators_evaluated = 1;
			data->People[next_individual].
			LS_performed[LS_2Opt_INDEX] = 0;
			data->People[next_individual].survival = 0;
		}
		indcpy(data, &(data->TheBestOne), &(data->People[0]));

		FreePopulation(&greedyPop);

#ifdef LOG
		/* Display partial results */
		log_on_file(data);
		dumpPopulation(data, solution_to_copy);
#endif
		for (int p = 0; p < parameters; p++) {
			data->stat.best_found_at_phase[0][p] =
			    data->TheBestOne.Score.measure[p];
		}
		/*
		if (data->BestMosInd.Score.measure[MOS] <=
		    data->MOS_Lower_Bound)
			return;
		*/

		/* Migliora gli individui generati con algoritmi costruttivi */

		if (data->LocSearchFreq > -2) {
			data->CurrentPhase = -1;
			for (next_individual = 0;
			     next_individual < solution_to_copy;
			     next_individual++) {
				ls_2opt_ind(data, next_individual);
				data->People[next_individual].survival = 0;
			}
			data->LB[TOTAL_COST] = 0;
			data->LB[VEHICLES_NUMBER] = 0;
			data->UB[TOTAL_COST] = INT_MAX;
			data->UB[VEHICLES_NUMBER] = data->CustomerSet.size();

			for (next_individual = 0;
			     next_individual < solution_to_copy;
			     next_individual++) {
				for (int p = 0; p < parameters; p++) {
					if (data->LB[p] >
					    data->People[next_individual].Score.
					    measure[p])
						data->LB[p] =
						    data->
						    People[next_individual].
						    Score.measure[p];
					if (data->People[next_individual].Score.
					    measure[p] > data->UB[p])
						data->UB[p] =
						    data->
						    People[next_individual].
						    Score.measure[p];
				}

			}
			update_k_lexicography(data);
			for (next_individual = 0;
			     next_individual < solution_to_copy;
			     next_individual++)
				fitness_function(data,
				                 &(data->
				                   People[next_individual].
				                   Score), data->k, data->LB,
				                 data->UB);
		}
		// riodina
		SortPopulation(data->People, solution_to_copy,
		               solution_to_copy);

		// TODO: aggiornare al caso corrente
		/*data->stat.best_mos_found_at_phase[1] =
		    data->BestMosInd.Score.measure[MOS];*/

#ifdef LOG
		/* Display partial results */
		log_on_file(data);
		if (data->LocSearchFreq > -2)
			dumpPopulation(data, solution_to_copy);
#endif
	}
#pragma endregion

#ifdef ORIG_SEQ
	/* Copia la sequenza originale... */
	onebody = data->People[next_individual];
	/* ...associando ai geni i valori [0 1 2 3 ... n] */
	for (g = 1; g < nGenes + 1; g++) {
		onebody.Sequence[g] = g;
	}
	for (int ls = 0; ls < LS_TYPES; ls++)
		onebody.LS_performed[ls] = 0;
	next_individual++;
#endif
	/* Allocazione degli individui random */
	for (; next_individual < nIndividuals; next_individual++) {
		onebody = &(data->People[next_individual]);
		/* Ad ogni individuo associa i geni [1 2 3 ... n] */
		for (g = 1; g < nGenes + 1; g++) {
			onebody->Sequence[g-1] = g;
		}
		/* Swappa a caso due geni */
		for (g = 1; g < nGenes; g++) {
			/* un numero a caso fra 0 e nGenes - g + 1 */
			x = rand() % (nGenes - g + 1);
			temp = onebody->Sequence[x];
			onebody->Sequence[x] = onebody->Sequence[nGenes - g];
			onebody->Sequence[nGenes - g] = temp;
		}
#ifdef _DEBUG
			for (int i = 0; i < data->NumGenes; i++) {
				assert(onebody->Sequence[i] != 0);
			}
#endif
		for (int ls = 0; ls < LS_TYPES; ls++)
			onebody->LS_performed[ls] = 0;
		onebody->performance_indicators_evaluated = 0;
		onebody->survival = 0;
	}
}

int stopEvolution(GeneticData * data)
{
#ifdef LOG
	log_on_file(data);
#endif

	/* check if the solution can be proved to be optimal */
	/*
	if (data->TheBestOne.Score.measure[TOTAL_COST] <= data->MOS_Lower_Bound) {
		data->stat.best_mos_found_at_phase[data->CurrentPhase + 2] =
		    data->BestMosInd.Score.measure[MOS];

		return AGA_EMOSLB;
	}
	*/

	/* check if the iterations (generations) limit has been exhausted */
	if (data->CurrentGeneration >= data->NumGenerations) {
		return AGA_EITLIM;
	}

	/* check if the time limit has been exhausted */
	/* if (data->timeLimit > 0
	   && difftime(time(NULL), data->beginTime) > data->timeLimit) {
	   return AGA_ETMLIM;
	   } */

	return 0;
}

inline void indcpy(GeneticData * data, Individual * to, const Individual * from)
{
	for (int i = 0; i < data->NumGenes; i++)
		to->Sequence[i] = from->Sequence[i];
	for (int j = 0; j < LS_TYPES; j++)
		to->LS_performed[j] = from->LS_performed[j];
	to->Score = from->Score;
	// gli open stacks? li copiamo o perderemmo solo tempo?
	to->performance_indicators_evaluated =
	    from->performance_indicators_evaluated;
	to->survival = from->survival;
}

/*
void mos_tos(GeneticData *data, Gene *sequence, int *mos, int *tos) {
	int *production, *stacks;
	int open, k, i, j;
	int n=data->NumGenes;
	int m=data->NumPieces;
	Pattern current;

	*mos = *tos = open = 0;
	production = (int *)calloc(m, sizeof(Piece));
	stacks = (int *)calloc(m, sizeof(Piece));

	// Setta la produzione al passo 0 e open e tos
	current = data->DM[sequence[0]];
	for (j=0; j<m; j++) {
		stacks[j] = production[j] = (current[j])?1:0;
		open += (current[j])?1:0;
	}
	*tos = *mos = open;

	// Chiude gli stacks da chiudere al passo 0
	for (j=0; j<m; j++) {
		if (production[j] == data->Cardinality[j] && data->Cardinality[j]) {
			stacks[j] = 0;
			open--;
			production[j] = 0;
		}
	}

	for (k=1; k<n; k++) {
		i = sequence[k];
		current = data->DM[i];

		for (j=0; j<m;j++) {
			if (current[j]) {
				if (!stacks[j])
					open++;
				(*tos)++;
				stacks[j] = 1;
				production[j]++;
			} else if (stacks[j])
				(*tos)++;
			if (open > *mos)
				*mos = open;
		}

		for (j=0; j<m;j++) {
			if (production[j] == data->Cardinality[j]) {
				stacks[j] = 0;
				open--;
				production[j] = 0;
			}
		}
	}

	free(production);
	free(stacks);
}
*/

//void print_step(int *stacks, int *production, int open, int tos, int m, int step, char *test_phase, int *Cardinality, int *current) {
//      printf ("Passo %d: %s\n", step, test_phase);
//      printf ("           current     [ ");
//      for (int j=0; j<m-1;j++)
//              printf("%d ", current[j]);
//      printf("%d ]\n", current[m-1]);
//      printf ("           stacks      [ ");
//      for (int j=0; j<m-1;j++)
//              printf("%d ", stacks[j]);
//      printf("%d ]\n", stacks[m-1]);
//      printf ("           production  [ ");
//      for (int j=0; j<m-1;j++)
//              printf("%d ", production[j]);
//      printf("%d ]\n", production[m-1]);
//      printf ("           cardinality [ ");
//      for (int j=0; j<m-1;j++)
//              printf("%d ", Cardinality[j]);
//      printf("%d ]\n", Cardinality[m-1]);
//
//      printf ("                 ( open, tos  ) = ( %d, %d )\n", open, tos);
//}
